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ABSTRACT 

We present a detonating failed deflagration model of Type la supernovae. In this model, the 
thermonuclear explosion of a massive white dwarf follows an off-center deflagration. We conduct a 
survey of asymmetric ignition configurations initiated at various distances from the stellar center. 
In all cases studied, we find that only a small amount of stellar fuel is consumed during deflagration 
phase, no explosion is obtained, and the released energy is mostly wasted on expanding the progenitor. 
Products of the failed deflagration quickly reach the stellar surface, polluting and strongly disturbing 
it. These disturbances eventually evolve into small and isolated shock-dominated regions which are 
rich in fuel. We consider these regions as seeds capable of forming self-sustained detonations that, 
ultimately, result in the thermonuclear supernova explosion. 

Preliminary nucleosynthesis results indicate the model supernova ejecta are typically composed of 
about 0.1 — 0.25 Mq of silicon group elements, 0.9 — 1.2 Mq of iron group elements, and are essentially 
carbon-free. The ejecta have a composite morphology, are chemically stratified, and display a modest 
amount of intrinsic asymmetry. The innermost layers are slightly egg-shaped with the axis ratio 
ft! 1.2 — 1.3 and dominated by the products of silicon burning. This central region is surrounded by a 
shell of silicon-group elements. The outermost layers of ejecta are highly inhomogeneous and contain 
products of incomplete oxygen burning with only small admixture of unburned stellar material. The 
explosion energies are « 1.3 — 1.5 x 10^^ erg. 

Subject headings: hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: 
general 



1. INTRODUCTION 

Almost half a century ago. lHovle fc Fowleil (|1960f ) pro- 
posed that some supernovae may originate from the de- 
generate remnants of stellar evolution. These objects are 
known as Type la supernovae (SN la) and are commonly 
believed to be the end points of the evol ution of interme- 
diate mass stars in close binary systems (jWhelan fc IbenI 
[T973h . 

The two most attractive theories of formation of 
Type la supernovae are the s ingle- degenerat e (SD; 
Whekn fc Ib^l [l973l: iN omotol [1981 IStarrfield et all 
2004t lYoon. Langer. fc Scheithauer 2005) and double- 



degenerate fDD: llben fc Tutukov. ,1984 iWebbinkI 1 19841) 
scenarios. The observational evidence necessary to dis- 
criminate which formation channel is preferred in nature 
remains indir e ct an d fragmentary (Branch c t al. 1 995t 
Livio fc Riesi l2003t iMannucci. Delia ValleT fc F'anagial 



20061 and references therein), in striking contrast to that 
of Type II supernovae. Evidence supporting the SD 
scenario has been collected only recently and in some 
cases requires more careful analysis. Some examples 
include observations of accreting massive white d warfs 
(|Lanz et all 120051 : ISuleimanov fc Ibragimovl l2003f) . evi- 
dence of hydrogen-rich material in the vicinity of the 
supernova, likely associated with a companion star 
hamuy-t-03, and the presence of the fast moving non- 
degenerate star inside a post-Type la supernova rem- 
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nant ()Ruiz-Lapuente et al.ll200^ . a possible companions 
of the supernova progenitor. In the case of the DD for- 
mation channel, several close binary white dwarf systems 
with total mass around Chandrasekha r mass have been 
identified during the last few years (|Napiwotzki et al.l 
120031 l2005l ). essentially confirming double degenerates as 
prospective progenitors of thermonuclear supernovae. 

Here we limit our considerations to a single-degenerate 
scenario in which the ignition of the degenerate mat- 
ter takes place in the core of a Chandrasekhar-mass 
white dwarf. Only a small difference is expected be- 
tween the outcomes of the core ignition of a mas- 
sive white dwarf formed through either the SD or DD 
channel. In the latter case, the explosion is perhaps 
born in the core of a massive rotating white dwarf 
(jPiersanti et al.l I2003D , a remnant of the final merger 
hase that does not result in an instantaneous explosion 
Guerrero. Garcia-Berro. & Isern 2004). 

In defining the initial conditions for multidi- 
mensional hydrodynamic models, we were guided 

in this study by res ults of t he recent analytic 

llGarcia-Se nz fc WoosleV '1995^; 'Wunsch fc WooslevI 
[2004!; 'Woosley, Wunsch, fc Kuhlcn 2001 and 
preliminar y numerical (|Hoflich fc Steinl 120021 : 
[Kuhlen. Wooslev. fc Glatzmaieil I2006D studies of 
conditions prevailing in the white dwarf's core just prior 
to the thermonuclear runaway. Our limited knowledge 
of that evolutionary phase grants us certain freedom 
in defining starting models. Both the number of the 
ignition points and the timing of the ignition are free 
parameters of the current models. 
Following the failure of carbon (jArnettj 119691 : 
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lArnett. Truran. fc Wooslevl 119711: lArnettI 11974) an d 
helium detonation models (jLivne k. ArnettI Il995f ). 
explosion modeling has focused on deflagration 
models and their derivatives. This preference has 
been firmly established after the apparent suc- 
cess of parametri zed one-dimensional pure defla- 
gration models (iNomoto. Sugimoto. fc Ned 119761 : 



iNomoto. Thielemann, fc Yokoil I1984D . They, along 
with t he later variant known as a de l ayed-detonation 
model ()Khokhlovl [l99l IWooslevI [T990t lArnett fc Li^ 
ll994bD brought the parametrized models into qualita- 
tive agreement with observations (jHoflich fc Khokhlovl 
[19961). 

Multi-dimensional studies initially included both 
2-dimension a l mode l s of d e flagration s ("Miillcr fc Arnetj 
[1981. [19861 : [Lh^ [l99l ILivne fc Arnet t 1995) and 
detonations ( Livne fc ArnettI Il995f ). During the 
last decade, sophisticated 3-dimensional models 
ha ve doni in ated t he scene. Central single-point 

fhokhlov '2000': 'Rcincckc, Hillebrandt. fc Niemeveii 
)2; Ga mezo et ah. 2003) as well as mult i-point igni- 
tion deflagrations 
20021 iRopke et al 



Garcia-Senz fc Bravol l2005h 



t;incckc, Hillebrandt, fc Nicmcver 
1 



2006.: 



.Schmidt fc Niemever, ,2006 : 
produce sublu- 



seem to 

minous events with highly mixed ejecta. Although 
models with a par ametrize d deflagration - to-det o nation 
transition (iGamez o. Khok hlov. fc OranI |2004 120051 : 
iGolombek fc Niemever. ,20051 ) appear to address both 
deficiencies, the mechanism behind the transition to 

detonation demands explanation; 

There exists e vidence (iHoflich et all l2002l 12004 
iKozma et al]|2005l but see also iBlinnikov et all (|2006f )) 
that the compositional structure of the ejecta ob- 
tained in multi-dimensional centrally-ignited deflagra- 
tions may not be compatible with observations. At 
the same time, spectroscopic and polarization obser- 
vations of several SNe la suggest the existence of 
strongly inhomog eneous outer ejecta layers rich in inter- 
mediate elements (IChugai. Chevalier, fc Lundavistll2004 
Wang et all l2004 iKasen fc Plewal 120051 : [Leonard et all 



20051 ) . These two apparently contradictory requirements 



can possibly be reconciled within a class of hybrid mod- 
els that combine deflagration and detonation within a 
single evolutionary sequence. Here we introduce a deto- 
nating failed deflagration (DFD) model, a generalization 
of ou r early exploratory study (jPlewa. Calder. fc LambI 
|2004 hereafter PCL), in which both essential elements 
are naturally present. In this model, the inhomogeneities 
in the outer ejecta layers result from the large scale per- 
turbation of the surface stellar layers induced by an off- 
center deflagration that fails to unbind the star. In the 
long term, that perturbation eventually leads to the for- 
mation of isolated shock-dominated regions that serve 
as ignition points for a detonation. The resulting event 
is luminous with a composite ejecta structure consist- 
ing of smoothly distributed detonation products in the 
central regions surrounded by inhomogeneous, turbulent- 
like outer layers composed of material partially burned 
in the deflagration. 

Our numerical investigations include certain simplifi- 
cations with the assumption of axisymmetry being po- 
tentially the most important one. This, however, will 
allow us to conduct a small parameter study exploring 
the dependence of the explosion parameters on the initial 



conditions. In turn, for the first time, we will be able to 
offer evidence that the simplifying assumptions regard- 
ing the geometry may not be the major deficiency of our 
models. 

2. METHODS 

We study the hydrodynamic ev olution of a massiy e 
white dwarf using the flash code (jFrvxefl et al.ll200d ). 
FLASH has been the subject of rather extensive verifi- 
cation a nd validation in both subsonic and supersonic 
regimes (| C alder et al.ll2002l: [Weirs et al.ll2005[ ). We used 
a customized version of the code based on the flash2.4 
release with specialized modules designed to model de- 
flagrations and detonations. We recorded the history of 
individual fluid elements with tracer particles for the pur- 
pose of future detailed nucleosynthesis studies. 

2.1. Reactive hydrodynamics 

We solved the time-dependent reactive Euler equations 
of self-gravitating flow in cylindrical geometry assum- 
ing axial symmetry. The non-reactive set of equations 
were extended by an advection-diffusion-reaction (ADR) 
equation describing the evolution of a deflagration front. 
The solution to the ADR equation w as obtaine d with 
the help of a flame capturing method (|Khokhlovl ll995). 
The FL ASH implementation has been the sub ject of verifl- 
cation (jVladimirova. Weirs, fc Rvzhikll20^6l ) with the re- 
sults of the application to turbulent fl ames closely rn atch- 
ing the original implementation ( Zhang et al.ll200^ . 

2.2. Flame model revision 

For the present application, several elements of the 
original ADR scheme were modified. In particular, care- 
ful analysis of the original three-stage FLASH burner 
(jCalder et al.ll20d^ PCL) revealed that it overestimated 
the amount of energy produced by burning the stellar 
mix to nuclear statistical equilibrium (NSE) at densi- 
ties typical of the stellar core. As a consequence, the 
nuclear ashes were both too hot and too rarefied, with 
buoyancy effects overestimated by a factor of w 3. This 
led to a much larger acceleration of the burned material, 
shorter evolutionary time scales (from ignition to bub- 
ble breakout), a lower consumption of stellar fuel and a 
correspondingly lower degree of pollution of the surface 
layers. The current scheme captures the energetics of the 
defiagrating material more closely. We also introduced a 
revised formula for the lam inar flame speed that bett er 
approximates the results of iTimmes fc Wooslevl ()1992[ ). 

2.2.1. Thermonuclear burning 

Following 'Khokhlovl (|199lL l2000f ). the evolution of a 
deflagrating stellar material can be considered as a se- 
quence of largely independent processes. Carbon burn- 
ing (stage I) is followed by relaxation toward nuclear sta- 
tistical quasi-equilibrium (NSQE) that produces silicon- 
group nuclei (stage II). Eventually, the matter relaxes 
toward nuclear statistical equilibrium (NSE) producing 
iron-group nuclei (stage HI). In the approximate burner, 
this is modeled through modifying the composition in all 
three stages. Stage I converts carbon into ^*Mg while 
^^Si is produced in stage II by the "burning" of ^^O and 
^''Mg. We also extended the original scheme by intro- 
ducing a "light nucleus" that accounts for the presence 
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Fig. 1. — Dependence of the energy release on the initial tem- 
perature obtained with torch47 nuclear network for a 50/50 C/O 
mixture over a fixed amount of time, tf = Tg. The highest 
amount of energy is released for Ti = 2 X 10® K and the lowest 
for Ti = 1.5 X 10® K. The energy release is essentially insensitive 
to the initial temperature for densities ^ 1 x 10^ g cm~^. 
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Fig. 2. — Temporal evolution of the energy release in self-heating 
torch47 nuclear network calculations for a 50/50 C/O mix at tem- 
perature Ti = 1.7 X 10® K. The bulk of the energy is produced 
within a single hydrodynamic time scale, Tg . Thick solid lines cor- 
respond to the energy release adopted in our supernova explosion 
calculations. 



of free alpha particles and protons in the very high den- 
sity and temperature regime (stage III, see below). 

The original flash approximate burning scheme as- 
sumed the energy release was a simple sum of the nom- 
inal differences between the binding energy of the ini- 
tial C/O mix and the burned material, Ai?b « 7.8 x 
10^^ erg g~^, independent of density. This is over 70% 
more than the energy release obtained by iKhokhlovl 
(|1983l Fig. 2, AEb « 4.5 x lO" erg g-^) for the den- 
sities of interest during the initial explosion stages, p « 
2 X 10^ g cm~^. Our revised sc heme uses resu lts ob- 
tained with the torch47 network ()Timmesll2001[ ) in the 
self-heating mode. The outcome of such calculations de- 
pends on the initial temperature of the fuel. Calculations 
also need to be conducted for long enough to guarantee 
a complete burn. 

In our first set of calculations, we advanced the net- 
work for a fixed amount of time, tf = Tg, where Tg = 
446p~^/^ s is the hydrodynamic ( free-fall) time scale 
(jFowler fc Hovldll964l : lArnettI 119961 ) . Calculations were 
started using several different initial temperatures of the 
mixture, Ti = 1.5 — 2.0 x 10^ K. Wc found a relatively 
weak dependence of the energy release on the initial tem- 
perature (Fig. [T|) and adopted Ti = 1.7 x 10^ K in sub- 
sequent calculations. 

Once the initial temperature was fixed, we turned our 
attention to studying the time-dependence of the energy 
release. We found that for p > 1 x 10^ g cm"'^ the bulk 
of the energy is released within At = Tg (Fig. [J). In the 
supernova calculations that follow, we used the energy 
release obtained for Ti = 1.7 x 10^ K and tf = Tg (thick 
solid line in Fig. [5]). The energy release was tabulated 
as a function of density with resolution 0.1 dex and lin- 
early interpolated. The flame extinction was modeled 
by smoothly decreasing the original energy release for 
densities < 1 x 10^ g cm~^. No energy was released for 
densities below 1 x 10^ g cm~^. 

In contrast to the original flash implementation of 
stage III in which ^^Ni was the sole product of burn- 
ing, here we introduced additional "light nuclei" com- 
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Fig. 3. — Composition of the "light nuclei" supplementing the 
NSE composition in stage III of the approximate burning scheme. 



posed of alpha particles and protons. This allowed us 
to better approximate the temperature, and in turn the 
buoyancy and dynamics of the Rayleigh- Taylor unstable 
burning front, especially at early times. As before, we 
used the torch47 nuclear network in an isochoric self- 
heating mode of burning to determine the composition 
of the light nuclei (Fig. [3]) . These results appeared insen- 
sitive to the particular choice of the initial temperature 
or the final time provided the network was evolved for at 
least Tg. 

2.2.2. Coupling to Hydrodynamics 

In addition to tracing the compositional evolu- 
tion of the nuclear ashes, the implementation of the 
Khokhlov's three-stage burner involves the coupling of 
the energy source term to hydrodynamics. In what 
follows we are primarily concerned with the early 
evolutionary stages when the physical time scales 
associated with nuclear burning are much shorter 
than a flame crossing time for computational zones 
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(iKhokhlov. Miiller. fc Holichl fl99l lArnett fc Livnd 
Il994bl : iReinecke. Hillebrandt. &: Niemeveill2002l ). When 
the evolution of the burning region is relatively slow and 
the expansion of the star is insignificant, as it is during 
the early evolution of centrally ignited deflagrations, a 
step-wise form of the energy production rate leads to 
large (about two orders of magnitude) rapid fluctuations 
in the energy deposition. These localized discrete 
"explosions" occur in partially-burned material where 
degeneracy is lifted and are a source of pressure waves 
and velocity fluctuations of order ~ 10 km/s. A similar 
problem in th e context of modeling deflagr a tions fronts 
was noted by IKhokhlov. Muller. fc HoflichI (|1993l ) who 
used, however, a completely different procedure to 
model the flame evolution. The origin of the problem, 
i.e. discrete representation of a thin unresolved burning 
front, is common to both studies. 

The impact of the numerical artifact just described can 
be limited by appropriately scaling the energy deposition 
in stage II (no such procedure is needed in stage I which 
always operates under strongly degenerate conditions 
and no energy is released in stage III). In the simula- 
tions presented here, the energy generation rate for stage 
II was kept at 1% of its nominal value throughout sim- 
ulations. Given the advection time scale « 10"'' s, this 
procedure affected the energy release only at very high 
densities for which the problem of spurious numerically- 
induced flow perturbations was originally identified. Al- 
ternatively, one can also consider less intrusive proce- 
dures in which the energy generation rate is limited only 
for times < 0.5 s. Imposing less strict limitations on the 
energy release was found to universally produce velocity 
perturbations of magnitude comparable to that of large 
scale flows. 

2.2.3. Flame speed calibration 

Several factors may cause the numerical flame speed 
propagation to differ fror n the prescribed laminar flame 
speed in one dimension. IVladimirova. Weirs, fc RvzhikI 
()2006[ ) have investigated in some detail the influence of 
numerical resolution, evolution with a superimposed con- 
stant background velocity mimicking hydrodynamic ad- 
vection (i.e. to verify Galilean invariance), and the im- 
pact of velocity gradients across the flame front (presum- 
ably caused by thermal expansion in real applications). 

In the supernova context, given the degenerate na- 
ture of the equation of state, thermal expansion is con- 
trolled by the energy release which chiefly depends on th e 
fuel density and composition (jTimmes fc Wooslevlfl99l) . 
Therefore whenever the flame energetics are modified, an 
appropriate calibration of the numerical parameters con- 
trolling the numerical flame speed needs to be done for 
the numerical flame speed to match the nominal (lami- 
nar) flame speed. 

The calibration procedure is relatively simple, but te- 
dious because the numerical flame speed depends on den- 
sity, composition, and numerical resolution. Given our 
focus is on the evolution of progenitors composed of a 
50/50 carbon/oxygen mix, our flame speed calibration 
was limited to that composition. We performed a num- 
ber of one-dimensional flame propagation simulations in 
Cartesian geometry. Models were obtained on grids with 
lengths between IL and 30i with L — 480 km, and us- 
ing between 256 and 1024 zones. The typical grid reso- 



lution was « 1 km. The results showed only weak de- 
pendence on the resolution. The flame propagation was, 
however, somewhat vulnerable to background flow fluc- 
tuations generated by the flame front motion that freely 
propagated and partially reflected off the boundaries. Al- 
though we used non-reflecting (zero gradient) boundary 
conditions, such reflections are understandable given the 
subsonic nature of the problem. In our analysis we used 
data from carefully selected long evolutionary sequences 
computed on the largest grids and at the highest possible 
resolution. 

Our calibration procedure was applied toward 
a slightly modified form ula originally proposed by 
iTimmes fc WooslevI (|1992l Eq. 43): 

si,TW = apg X (^imi - e"(^) ^ 

with a = 35.46538 x 10^ b = 1.110603, Sq = 2.6132427 x 
10"^ and ag = 2.9538546 x 10' ^. This improved formula 
reproduces iTimmes fc WooslevI (|l992) results for a C/0 
50/50 mix to within 15% for densities between 1 x 10^ and 
1 X 10^ g cm"-^. Finally, this has been adjusted to account 
for thermal expansion due to the torch47 energetics: 

SLFLASH = si^TW X max(0.9, min(1.3, p(logio p))), 
where 

p{x) — Cq + ClX + C2X^ -f CzX^ -I- c^x^ , 

is polynomial density-dependent correction factor based 
on our calibration calculations and cq = 413.6563, 
ci = -194.3208, C2 = 34.06912, cg = -2.633218, and 
C4 = 7.5665459 x 10~^. The final formula reproduces 
ITimmes fc WooslevI (|1992l Table 3) resuhs to within 5% 
for densities in the range 1 X 10'^ and 4 x 10^ g cm . 

2.2.4. Hybrid burning scheme 

Our initial investigations into late evolutionary stages 
of fizzle off-center defiagrations (PCL) provided strong 
evidence that the conditions inside the confluence re- 
gion at the stellar surface are appropriate for creating 
a detonation. We observed that both densities and tem- 
peratures were high enough for the burning time scale 
to become shorter than hydrodynamic time scale. We 
also anticipated that the shocked region would remain 
confined for an extended period of time sufficient to 
develop a self-sustained detonation. We were unable, 
however, to study that process in detail at that time. 
The burning module could not reliably discriminate be- 
tween shocked fuel (a legitimate detonation site) and 
compressed partially burned matter. Here we introduce 
a hybrid burning scheme to allow for the simultaneous 
presence of deflagrations and distributed nuclear burn- 
ing in the simulation. Defiagrations are modeled us- 
ing the ADR flame-capturing schem e (|Khokhlovl 119951 : 
IVladimirova. Weirs, fc RvzhikI [200^ with modiflcations 
as outlined above. The distributed burning is calcu- 
lated using the flash aproxl3t 13-isotopc alpha net- 
work (G. Jordan 2005, private communication) This net- 
work is an extension of the original flash approxlS 
network and includes temperature coupling for increased 
stability of calculations in the NSE regime Miillcr (1986). 

The first step in our our hybrid burning procedure is 
to identify shocked zones. This is done using a multi- 
dimensional shock detection module adopted from the 
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SPPM code (| Anderson fc WoodwardlfTQQSD with pressure 
jump across the shock, Ap/p > 0.5. If the tempera- 
ture in such zones is high enough for nuclear burning 
{Tnuc,min = 8 X 10^ K), a progress variable used by 
the ADR flame module is reinitialized and the numer- 
ical flame speed is set to zero. The flame module is 
subsequently activated only if the flame speed is greater 
than zero and the material is pure fuel (as indicated by, 
for example, a small abundance of nuclei not participat- 
ing in the simplified three-stage burner). Otherwise, the 
nuclear network is called but only if the following con- 
ditions are fulfilled simultaneously: (1) both the density 
and temperature are in the range appropriate for the net- 
work calculations, (2) the zone is outside a shock front, 
(3) the defiagration module was not used, and (4) the 
ADR progress variable is small. The first condition limits 
network calculations to high-density (p > 1 x 10^ g cm^'^) 
high-temperature (T > 8x 10^ K) regimes. Condition (2) 
follow s from the conclusions of lFrvxell. Miiller. fc Arnefrg 
([1989') who found that the correct speed of detonation 
waves can be obtained only in models with burning dis- 
allowed inside unrelaxed numerical shock profiles. Con- 
ditions (3) ensures only one physical burning process, 
either distributed burning following shock heating or de- 
flagration, operates at a time. Finally, condition (4) pre- 
vents burning from occurring in regions preheated by 
an extended diffusive tail of the ADR flame capturing 
scheme and at the same time allows for burning to con- 
tinue after the detonation wave sweeps through partially 
deflagrated material. By making the progress variable 
threshold for use of nuclear network slowly increasing 
with time, we also prevented spurious (distributed) burn- 
ing in the dense central regions perturbed by ascending 
deflagrating material. 

The above scheme captures the essential behavior of 
both deflagrations and wave-induced distributed burn- 
ing. It allows for the transition to detonation in shock- 
heated regions away from deflagrating material and en- 
ables burning in shocked, partially deflagrated material. 
The above selection rules were developed and improved 
based upon experience accumulated in the course of sev- 
eral numerical experiments and, as such, offers a prac- 
tical rather than ideal recipe for modeling deflagrations 
and detonations in the same physical setting including a 
possible deflagration to detonation transition. 

Once the detonation wave is launched, a combina- 
tion of a PPM hydrodynamic solver, a multidimensional 
shock detection algorithm, and a aproxl3t nuclear burn- 
ing module was used to advance the evolution in time. 
The evolution of the detonation wave on a large scale is 
expected to be captured correctly given that the thick- 
ness o f the wave is much smaller than the stellar radius 
()Fallel[200ft) . We did not find it necessary to rescale nu- 
clear reaction rates in order to o btain the correct propa - 
gation of the detonation speeds (lArnett fc Livnd[l994b( ): 
excluding the unrelaxed shock profile from burning ap- 
peared sufficient to yield physical solutions. 

2.3. Some comments on numerical model limitations 

One of primary motivations behind development of a 
hybrid burning scheme was desire to mitigate a risk of 
spurious detonation ignitions. The existence of such spu- 
rious ignitions is a well-documented fact in a strophysi- 
cal literature (|Frvxell. Miiller. fc Arnettlll989f) . Here we 



only briefly discuss the most typical causes for spurious 
ignitions and possible ways to prevent them from pollut- 
ing hydrodynamical models. 

Perhaps the most common cause for spurious deto- 
nation ignitions is a numerically-induced mixing of hot 
ashes and cold fuel. By advecting a material inter- 
face (contact discontinuity) sep arating ashes from fuel, 
iFrvxell. Miiller. fc ArnefrB (|1989D demonstrated that such 
mixing may result in artiflcial preheating of fuel, its ig- 
nition, and ultimately formation of a combustion wave. 
Related to this is a problem of species conservation by 
nonlinear Eulerian advection schemes. Under certain 
conditions numerical modification of fuel (or partially 
burned material) composition may change burning en- 
ergetics. Either extinction or enhancement of burning 
may follow. To our knowledge, neither numerical species 
diffusion nor species non-conservation can be completely 
eliminated from Eulerian simulations of realistic nonlin- 
ear systems. Our hybrid burning scheme attempts to 
limit possible infiuence of both effects by constraining 
burning to regions occupied by pure fuel. 

Transition to detonation can also follow an artifi- 
cial boosting of an acoustic perturbation. Such per- 
turbations lead to local variations in temperature and 
under degenerate conditions temperature is a sensitive 
function of (relatively small thermal) pressure. Given 
strong dependence of nuclear reactions on temperature, 
it is conceivable that even small but sustained heat- 
ing may strengthen acoustic waves and eventually cause 
spurious transition to detonations. Typically, however, 
small acoustic fiuctuations suffer from strong damping 
by numerical diffusivity and such spurious transitions 
to detonations are likely to occur only if nuclear burn- 
ing is allowed inside numerical (unrelaxed) shock pro- 
files. The cure for this problem is to eliminate burning 
in reg ions occupied by shocks (Fryxell, Miiller, fc Arnetfl 
Il989| ) and, as we mentioned above, such a filter is em- 
ployed in our calculations. 

Finally, application of nuclear burning source term in 
our calculations is limited to regions with sufficiently 
high densities and temperatures. That is, we wish to con- 
sider a feedback from nuclear burning only if the nuclear 
timescale is short enough to infiuence hydrodynamics. 
This approach not only saves computational time, but 
more importantly prevents the nuclear network from be- 
ing fed with input data representing low-density regions 
where evolution has a highly transient character and is 
not correctly captured in our calculations. 

To summarize, numerical computations of reactive 
flows pose extreme challenges and require very careful 
treatment. We attempted to address several known and 
some newly emerged problems related to coupling reac- 
tive source terms to hydrodynamics in great detail. The 
impact of some of these problems could only be limited, 
but not completely eliminated. For example, in our mod- 
els no nuclear burning is allowed inside numerical shock 
proflles. However, our particular choice of parameters 
deflning numerical shock profile may not be adequate 
in all situations affecting evolution of shocks and acous- 
tics in unwanted manner. Poor numerical resolution only 
adds to the algorithmic inefficiency further limiting pre- 
dictive abilities of our models. 

Clearly, successfully resolving technical problems of 
our computations is of high priority and such aspects 
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should be remembered when interpreting our results. At 
the same time, however, one should also keep in mind 
that our model involves several simplifying assumptions 
(i.e., we consider a non-rotating, non-convective, non- 
magnetized chemically homogeneous progenitor), and as 
such it is a blend of approximate numerics and unique 
choice of parameters defining physical scenario. 

2.4. SN la explosion code verification: central 
deflagrations 

Given that both the basic hydrodynamic module as 
well as the FLASH implementation of the flame captur- 
ing scheme have been extensively verified and, albeit to 
a lini i ted extent, a l so val i dated in the past (jCalder et al. 
2002t iWeirs et all l2005t iVladimiro va. Weirs, fc Rvzhik 
our verification is solely limited to code-to- 
code comparison in the context of thermonuclear super- 
nova explosion modeling. Although code-to-code com- 
parison is widely popular among computational physi- 
cists, the usefulness of this appr oach is hotly debated 
(jTrucano. Pilch, fc Oberkampfll2003i) . For one, even per- 
fect agreement between the simulation results of two 
codes does not offer proof of their correctness. More- 
over, the scope of such an exercise is usually limited by 
the specific capabilities of the code, the availability of 
the results, the completeness of documentation, and the 
relevance of the performed test to the actual problem at 
hand, to name a few. Here we accept the obvious de- 
ficiencies of the code-to-code comparison approach and 
use highly relevant and state-of-the-art calculations as a 
benchmark. Again, due to the novelty of our ultimate ap- 
plication, no data are available for comparison, though 
they will hopefully become available through indepen- 
dent calculations. Ultimately, the model will be validated 
using accumulated observational evidence. Methodology 
and methods required in such assessments are presented 
in a companion pape r (|Kasen fc Plewal 120061 see also 
iBlinnikov et"al] (|2006[ )). 

2.4.1. Simulation setup 

For our comparison study we selected a family of cen- 
trally ignited deflagration models obtained by the G arch- 
ing group (iReinecke. Hillebrandt. fc Niemeveil l2002t 
iRdpke fc Hillebrandti l2005b[ and references therein) . 
These studies are very well-documented and their results 
compare favorably to those obtained hy^ other group s 
(jGamezo et al.ll20'03HGamezo. Khoklflov.' fc Oranll200l . 
When comparing results, we considered the overall 
morphology of the explosion (flame front structure), 
the energetics and the approximate ejecta composi- 
tion of t he 2-dim ensi onal explosion calcu l ations re- 
ported bylReinecke. Hillebr andt. fc Niemeveil (|2002f ) and 
iRopke fc Hillebrandti (,2005^. Our choice of 2-D Garch- 
ing models is natural given our calculations to follow also 
assume axial symmetry. 

In a benchmark study we used flash2.4 customized 
for the thermonuclear supernova explosion problem. We 
used a PPM module for real gas inviscid hydrodynamics 
and the Helmholtz equation of state required by the de- 
generate conditions encountered in the white dwarf inte- 
rior. All calculations were done with Courant factor 0.6. 
This choice of time step limiter may appear somewhat 
conservative, but allows for better coupling between dif- 
ferent physical processes. 



Gas self-gravity was accounted for by solving the Pois- 
son equation through multipole expansion. We found 
that linear momentum is poorly conserved in explosion 
calculations when the expansion series is truncated too 
early, especially when the explosion displays significant 
asymmetries. In test calculations done assuming a spher- 
ical potential, the bulk of the stellar material displayed 
motion of a few hundred km/s after only a few seconds 
of evolution. The momentum conservation gradually im- 
proved as additional higher order terms were included in 
the expansion. In what follows we used 10 multipole mo- 
ments and found excellent momentum conservation for 
all initial conditions considered. 

We used a 2-dimensional cylindrical grid (r, z) and as- 
sumed axial symmetry. This implied imposing a reflect- 
ing boundary condition at rmin = 0. We applied outflow 
conditions at the remaining boundaries. In our verifi- 
cation calculations the computational domain covered a 
rectangular region with rmax = Zmax = -Zmin = 16, 384 
km. We used the adaptive capabilities of the flash code 
to create several levels of refinement up to a maximum 
resolution of 8 km. We do not expect the dynamical 
evolution of low density gas or at large distances from 
the stellar center to play an important role in explo- 
sion simulations and therefore adaptive refinement was 
used only for radii < 4, 000 km and if the gas density 
> 1 X 10* g cm"'^. Self-gravity calculations, on the other 
hand, require good resolution of dense regions and the 
grid refinement was forced to the highest allowed reso- 
lution whenever the gas density > 3 x 10^ g cm~'^. In 
addition, the innermost 2, 500 km of the star have al- 
ways been resolved with the finest zones. In regions 
where adaptivity was allowed, AMR patches were cre- 
ated dynamically if the local density contrast exceeded 
0.50 or the total velocity changed by more than 20%. 
Furthermore, we ensured the flame front was always re- 
flned to the flnest level by forcing reflnement whenever 
the fractional change of the flame progress variable ex- 
ceeded 0.02. 

2.4.2. Initial model 

The supernova progenitor was an isothermal, T — 
5 X 10^ K, white dwarf composed of equal mass frac- 
tions of carbon and oxygen. With a central density of 
2 X 10^ g cm~'^, the progenitor had a radius « 2100 km, 
total mass « 1.36 Mq, and total energy « -4.92 x 10^° 
erg. The progenitor was surrounded by a low den- 
sity {pamb = 1 X 10~^ g cm~^) and low temperature 
{Tamb = 3 X 10*" K) medium composed of pure helium. 
The stellar material and low density ambient medium 
were initially marked with a passively advected tracer 
that was set to 1 and in those two regions, respectively. 
Subsequently, gas gravitational accelerations were multi- 
plied by the tracer value in the course of the evolution, 
allowing us to prevent the ambient medium from collaps- 
ing onto the central object and limiting mass diffusion at 
the stellar surface. 

The original progenitor model was constructed using 
a numerical discretization different from the one used 
in the hydrodynamic simulations and assuming a sim- 
plified equation of state. For that reason, the perfect 
hydrostatic equilibrium of the original model was de- 
stroyed as soon as it was interpolated onto the simulation 
mesh. Even though the mismatch between the two com- 
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time [s] 

Fig. 4. — Evolution of the root mean square velocity in mod- 
els without burning obtained with resolution 2 (thick solid), 4 
(medium solid), 8 (thin solid), 16 (dotted), and 32 km (dashed). 
The initial models were perturbed by adding random velocities 
with amplitude 200 km s~^ inside the stellar core. Note that the 
velocity gently decays in models calculated with a resolution of 
16 km or better. 

putational environments was relatively small, very strong 
acoustic oscillations quickly developed making such a 
model unsuitable for further investigations. Moreover, 
the oscillations did not decay with time, presumably due 
to both low dissipation of the hydrodynamic scheme and 
the strong degeneracy of the medium. 

We constructed a stable progenitor model using a m od- 
ified variant of the damping method of lArnet^ ()1994[ ) in a 
1-D FLASH simulation in spherical geometry. Our proce- 
dure combined a very slow diffusion of velocity together 
with a partial rather than complete (as in the original 
method) removal of excess momentum after each time 
step. We found that the complete removal of momentum 
prevented the model from reaching equilibrium. Damp- 
ing process usually lasted several thousands of hydro- 
dynamic steps. We examined the stability of the relaxed 
model by computing a sequence of hydrodynamic models 
without burning. Random velocity perturbations with 
amplitude of 200 km s~^ (Ma « 0.02) were added to 
the inner core region of radius 400 km. We observed 
the decay of the root mean square velocity with time, as 
expected in a stable model, provided the resolution was 
16 km or better (Fig. We also computed the evolu- 
tion of an "effective" stellar radius corresponding to the 
volume occupied by gas with density > 1 x 10'' g cm~'^ 
representing the bulk of the stellar matter. The results 
are shown in Fig. O As one can see, the model stellar 
radius shows significant evolution in computations with 
resolution 16 km or worse; only a very modest (« 0.5%) 
decrease of radius was observed in the better resolved 
models, and had no consequence for structure of the stel- 
lar core (i.e. possible temperature increase). 

Unfortunately, close examination of the velocity field 
also revealed that while the overall stability of models 
improves with increasing resolution, small velocity per- 
turbations not only do not decay but are amplified near 
the symmetry axis. This effect was especially strong in 
2 km resolution model where the velocity near the axis 
rapidly increased from the initially imposed 200 km s^^ 
to 800 km s^^ during the first 0.1 s of the evolution 
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Fig. 5. — Evolution of the stellar radius in models without burn- 
ing obtained with resolution 2 (thick solid), 4 (medium solid), 8 
(thin solid), 16 (dotted), and 32 km (dashed). The initial mod- 
els were perturbed by adding random velocities with amplitude 
200 km s~^ inside the stellar core. Models computed with resolu- 
tion no worse than 8 km are stable, displaying only a very small 
degree of radius change. The model using 32 km resolution shows 
very strong expansion, making it unsuitable for long-term hydro- 
dynamic evolution studies. 
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Fig. 6. — Total velocity profiles at t = 0.1 s near the symmetry 
axis (r = cm) in models without burning computed with 2 km 
(thick line) and 8 km (thin line) resolution. The velocity ampli- 
tude rapidly and significantly increases in the 2 km model, while 
it slowly decreases in the 8 km model. 



(shown with the thin solid line in Fig. (5]). No signifi- 
cant increase in the magnitude of the spurious velocities 
was observed at later times, but the affected region ex- 
panded from the initial 1 to 3 zones by t — Is. We 
observed very similar behavior in the 4 km resolution 
model, although the velocities near the symmetry axis 
were ultimately somewhat smaller (« 600 km s~^). In 
contrast, the velocities smoothly decay from their initial 
values in the whole perturbed region in the 8 km resolu- 
tion model (shown with the thick solid line in Fig. [H]). A 
similar slow decrease of the velocities near the symmetry 
axis was also observed in the 16 km model. 

After a stabilized progenitor model was interpolated 
onto the simulation mesh, a deflagration front was ini- 
tialized around a small region at the stellar center. The 
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Centrally Ignited Benchmark Deflagration Models 
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flame front was located at 



Rf = R,n X 



1 



am cos 



Um tan 



where Rm is the unperturbed radius of the burned re- 
gion (flame radius), a„j is the flame radius perturbation 
amplitude, and n„ controls the flame radius perturba- 
tion wavelength. Random velocity perturbation were 
added to the inner core region following the procedure 
described above. By introducing a small stochastic com- 
ponent ("cosmic variance") into the problem, we were 
able to examine the robustness of our results. In par- 
ticular, it is essential to verify that the observables (i.e. 
the explosion energy) do not show strong dependence on 
small perturbations in the initial conditions, in accord 
with the observed intrinsic homogeneity of SN la. 

2.4.3. Centrally ignited benchmark models 

We performed a comprehensive survey of two- 
dimensional centrally ignited deflagration models. The 
database contains 33 models. All models were evolved 
until t = 2.5 s when burning essentially ceased. For each 
model we varied the flame perturbation wavelength and 
random velocity perturbation pattern (through the seed 
perturbation). We explored the sensitivity of the results 
to mesh resolution by resolving the central stellar region 
of radius Rc t Axc for times < tc- Resolution coarsening 
in this central region was done in equal intervals of time. 
For times > tc, the resolution was equal to a default value 
of 8 km. The initial flame radius and perturbation am- 
plitude were in all cases fixed at R„i — 100 km and 0.1, 
respectively. Table[T]presents a complete database of our 



benchmark models.^ The model identification tag name 
is constructed as a string n — d-r — t — ph where n — de- 
notes the perturbation pattern, d- describes the maximal 
grid resolution in the core region, r — denotes a radius 
of the core region inside which enhanced resolution was 
used, t — denotes the time up to which enhanced flame 
resolution was allowed, p distinguishes between differ- 
ent perturbation patterns, and finally h denotes the de- 
sired thickness of the flame front (in grid ce lls) in the 
ADR flame capturing module (see Appendix in 'Khokhlovl 
11995). The default value of the numerical flame thick- 
ness was 4; we obtained two models, nlld2r40t30b3 and 
nlld2r40t30b6, in which the nominal flame thickness was 
varied by —25 and -1-50%, respectively The last column 
in Table [T] gives the total energy of the model (sum of ki- 
netic, internal, and potential energies) at the final time, 
t = 2.5 s. 

Evolution m reference model — We use model 
nlld2rl0tl0a to illustrate major characteristics of 
a centrally ignited defiagration in our benchmark 
configuration (Fig. [7|. The density in this model at the 
initial time is shown in Fig. [7]Ja). The flame front has 
the shape of the initial regular n = 11 perturbation. 
At a resolution of 2 km, comparable to that of the 
highest resolution G arching group models at early times 
(|Ropke et al.ll2006l . Fig. 3), the flame region is resolved 
into about 50 zones in radius. This initial configuration 
undergoes a dramatic evolution during the next second 
(Fig. El^b)). After that time only 3 prominent bubbles 
are clearly identifiable and some parts of the fiame 
begins forming disconnected regions (e.g. the region 
located near (r, z) w (650,-550) km). The fiame leaves 
behind a significant amount of unburned material and 
a few isolated pockets of fuel can also be identified 
inside rising bubbles. The outermost parts of a highly 
convoluted flame front have reached ^ 1,000 km in 
radius. The star has expanded by « 20% and the typical 
expansion velocities are ~ 2, 500 km s^^. After another 
second (t = 2 s; Fig. EI^c)), the expansion of the outer 
stellar layers becomes slightly nonuniform due to the 
uneven acceleration caused by individual flame bubbles. 
The expansion velocity exceeds 10, 000 km near the 
stellar surface. At this time one can still identify 3 large 
bubbles, but these are now more developed and occupy 
a much larger volume fraction. Their morphology does 
not change much at still later times {t = 2.5 s, the final 
time; Fig. [Tjd)) as nuclear burning essentially quenches 
and the ejecta expansion becomes progressively radial. 

Sensitivity to small perturbations — Given the highly non- 
linear character of the Rayleigh- Taylor unstable defla- 
grating bubbles it is natural to ask how robust are the 
predictions offered by individual models. We studied this 
question by creating a sequence of models for a given set 
of primary model parameters (initial flame configuration 
and numerical resolution). The initial conditions for each 
member of a given sequence differed only by the pattern 
of small amplitude stochastic velocity added to the ini- 
tially static progenitor model. Our metric for compari- 
son is to examine the final ejecta morphology and integral 



^ The database is available on-line 
flash . uchicago . edu/tomek/Results/DFD/central/. 
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Fig. 7. — Hydrodynamic evolution of the centrally ignited benchmark supernova model nlld2rl0tl0a. Panels (a)-(d) show the density in 
logarithmic scale together with the outline of the flame front (contour line corresponding to progress variable value of 0.5). Notice the scale 
change between panels, (a) The initial conditions at t = s. (b) t = 1.0 s. The flame front is highly convoluted; the star remains spherical 
but expanded by 20%; the bulk expansion velocity is Ri 2,500 km s~^. {c)f = 2.0s. The flame front is rich in structure with some 
pockets of unburned material; a large amount of unburned material can be found near the center; individual flame bubbles create smooth 
large scale impressions on the stellar surface; the expansion velocity near the outer stellar edge are just slightly over 10,000 km s-l. (d) 
t = 2.5 s. The flame is extinct and the structure of the flame front becomes frozen. Evolution toward homologous expansion begins. 



model characteristics (temporal evolution of total energy, 
burning rate, flame surface area). 

As an example we compare select models from a single 
sequence. Figure [S] shows the ejecta temperature dis- 
tribution at the final time in the sequence nlld2rl0tl5. 
The isocontour of gas density of 1 x 10'' g cm~'^ is shown 
with the black line and can be identified with the stellar 
surface. Several comments can be made following inspec- 
tion of Fig. M 

It is encouraging to notice that the models show no 
axis-related bias, a numerical artifact that frequently 
pollutes axisymmetric hydrodynamic calculations. In 
particular, the structures developed near the symmetry 
axis (r = cm) in models nlld2rl0tl5a (Fig[5]Ja)) and 
nlld2rl0tl5b (Fig[Hb)) are markedly different. There 
is also no visible difference in the amount and quality of 
the structure developing in the regions above (z > cm) 
and below (2; < cm) the equator. Some models do 
develop structures near the equator (see, e.g., Fig[8ljd)), 
but some others do not (see, e.g., Fig[5]^c)). This allows 
us to conclude that possible defects due to the geome- 
try representation do not affect our calculations in any 
significant way. 

In all models considered here large scale structures 
(bubbles) several tens degrees in size dominate in the 
outer parts of the ejecta. In some models perhaps no 
more than 2 (Fig[H]Jc)) while in some others perhaps as 
many as 3 (Fig[5]^d)) such large and distinct structures 
survive turbulent burning. These bubbles push ahead 
unburned material causing relatively mild deformation 
of the ejecta outer layers as indicated by shape of the 
density isocontour. 

As discussed below, the explosion energies also ap- 
pear sensitive to small perturbations. In the case of 
the nlld2rl0tl5 models, our limited sample shows the 
total variation in explosion energy « 0.1 foe (1 foe = 
1 X 10^^ erg) or w 10% in the released energy (see 
Fig. IDJc)). This may indicate that small, naturally oc- 
curring variations in the internal structure of progenitors 
(expected to arise from the convective flows developing 
in their cores prior to runaway) may contribute to the 
intrinsic diversity of SN la. Addressing this interesting 



possibility requires more careful study, preferably using 
realistic multi-dimensional progenitor models. 

Sensitivity to numerical resolution — From the mod- 
eler's point of view, several simulation parameters may 
potentially affect the results and therefore need to 
be controlled. Given the high degree of complex- 
ity and highly non-linear character of our applica- 
tion, it is natural to expect that the model results 
will depend upon the numerical resolution. Conver- 
gence to the true 3-dimensional solution is not ex- 
pected to occur in two dimensions due to, for ex- 
ample, differences between the physics of tw o- versus 
three-dimensional supernova turbulence (iKhokhlov 19951: 
ISchmidt. Hillebrandt. &: Nie mcvcr 2006) and Rayleigh - 
Tavlor instabilitv p<ane et al.. .2000; .Chertkovl ]2003h . 
Nevertheless, it is still important to examine the sensitiv- 
ity of our axisymmetric model predictions to numerical 
resolution. 

In adaptive mesh simulations, the computational mesh 
is not necessarily a well-defined entity as the numeri- 
cal discretization depends on the solution and is usually 
highly variable both in time and space. The simple pro- 
cedure of doubling the grid resolution does not have its 
usual interpretation. Unlike uniform grid simulations, 
adaptive mesh refinement computations admit additional 
error into the solution by not resolving smooth or other- 
wise dynamically insignificant parts of the flow field. In 
addition, some errors, such as flow perturbations arising 
at the fine-coarse mesh interfaces, are uni que to AMR 
discretization and not easy to charact erize (|QuirkllT991l : 
IWeirs et al.1l2005l: iPantano et al.ll2005[ ). 

As our earlier investigations have demonstrated, the 
morphology on small scales of the exploding models ap- 
pears very sensitive to slight perturbations in the initial 
conditions, and might be useful only for making qualita- 
tive statements. A more quantitative comparison of the 
different models can be made using integral quantities. 
For example, variations in the final explosion energy are 
of great interest from the observational point of view, 
and several possible natural causes for such variations 
have been proposed (i.e. differences in the chemical com- 
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Fig. 8. — Select centrally ignited benchmark supernova models from sequence nlld2rl0tl5. Panels (a)-(d) show the temperature in 
logarithmic scale at the final time (t = 2.5 s) in models nlld2rl0tl5a, nlld2rl0tl5b, nlld2rl0tl5c, and nlld2rl0tl5d, respectively. The 
density isocontour at p = 1 X 10'* g cm""^ is shown as the black solid line. Although the details of ejecta morphology vary strongly between 
models, all models produce large bubbles of burned material gently deforming the outer stellar layers. The presence of a symmetry axis in 
the simulation does not appear to bias the calculations and no asymmetry between the two hemispheres is observed. 



position and/or the rotation of the progenitor). 

Here we discuss the role of mesh adaption on the evo- 
lution of the total energies, flame surface areas and the 
burning rates. Since as we mentioned before the evolu- 
tion is rather sensitive to small perturbations, we chose 
to compare different families of models rather than in- 
dividual family members. Figure [H] shows the evolu- 
tion of the integral quantities in three families of the 
nlld2 sequence of benchmark models: nlld2r05tl0 (left 
panel), nlld2rl0tl0 (middle panel), and nlld2rl0tl5 
(right panel). The individual members of each family 
were obtained using slightly different random velocity 
perturbations. Compared to the nlld2rlGtlG models, 
in models nlld2r05tl0 the innermost region of enhanced 
resolution (Ax — 2 km) was limited to 500 km in radius 
(as compared to 1 , 000 km) . The resolution in this region 
was decreased by a factor of 2 at 0.5 s and by another 
factor of 2 at t = 1.0 s. In the nlld2rl0tl5 family of 
models, this innermost region was derefined in two simi- 
lar stages with the resolution ultimately decreased to its 
nominal value at t = 1.5 s. 

Several observations can be made following analysis of 
the the nlld2 sequence. 

All nlld2 models produce explosions. The explosion 
energies are rather low, « 0.4 — 0.5 foe. In all cases 
the nuclear burning is most intense around « 1.2 s after 
the ignition. That period of rapid energy release lasts 
for about 0.5 s with the densities dropping rapidly to 
« 1 — 2 X 10'' g cm~^ hy t = 1.5 s. This is followed by 
a period of low, almost constant energy release during 
which in only exceptional cases a slight increase in en- 
ergy generation was observed. At that time (t sa 1.7 s), 
the burning takes place at densities < 1 x lO'' g cm^'^ and 
changes its character from active turbulent burning asso- 
ciated with vigorous creation of flame surface to a much 
milder, distributed mode of burning. By t = 2 s the 
densities drop to ~ 1 x 10^ g cm~'^ and flame quenching 
results in a steady decrease in energy generation. 

Models obtained with higher resolution produce more 
energetic explosions. The typical explosion energies vary 
from w 0.45 foe for the least resolved subsequence r05tlO 
to « 0.55 foe for the best resolved subsequence rl0tl5. 
Higher resolution also appears to make evolutionary tra- 



jectories more similar at early times (the curves run more 
closely in subsequence models rlOtlO than in r05tl0) but 
result in increasing diversity at late times (there are rela- 
tively large variations in energy generation rates around 
i = 1.5 s in subsequence rl0tl5). 

We found that not only the morphology, but also the 
integral quantities are sensitive to small perturbations 
in the progenitor. For example, the dispersion of explo- 
sion energies is about 0.1 foe, even in relatively well- 
resolved simulations (e.g. subsequence rl0tl5). This 
may indicate that some of the observed diversity of 
supernovae might be produced by the nonlinear re- 
sponse of the explosion process to small variations in 
the initial conditions. Such variations from one pro- 
genitor to another are expected to exist in nature 
especially given the convective (turbulent) flow con- 
ditions prevailing in th e stellar cores prior to ther- 
monuclear runaway (see 'Wooslev"l990'; 'Hofli ch fc Steid 
2002; Kuhlcn, Wooslcy, & Glatzmaicr 2006, and refer- 
ences therein). The contribution of such a purely 
stochastic component to the explosion process clearly de- 
serves more careful study. 

2.4.4. Comparison against Garching group models 

We found good agreement between the main charac- 
teristics of our centrally ignited model explosions and 
the results of equivalent axisymmetric calculations pre- 
sented by the Garching group. The overall evolu- 
tion of the energy generation rate in the nlld2 model 
sequence is similar to that in m odel c3_2d_256 by 
(jReinecke. Hillebrandt. &: Niemev cr 200^ Fig. 3) and 
to that obtained earlier bv (iNiemeveii Il995l ; see 
iReinecke. Hillebrandt. fc Niemeveil (|1999D ). The energy 



generation in our model explosions displays a pronounced 
maximum reaching between w 1.1 — 1.2 x 10^^ erg s~^ 
in models nlld2r05tl0 and nlld2rl0tl0. This com- 

J ares v ery favorably to the result reported b y iNiemeveil 
19951 ). IReinecke. H illebran dt. fc Niemeved (|2002D . and 
more recently by Ropkg (j2005l ). The latter two studies re- 
ported peak energy generation rates ~ 1.2 x 10^^ erg s~^. 
The rates obtained in models nlld2rl0tl5 are higher 
by about 50%. S uch higher rates were reported by 
(|Ropke et al.l l2006f ) for some of their centrally ignited 




Fig. 9. — Evolution of the total (explosion) energies and burning rates (solid and dashed lines, respectively) in three families of centrally 
ignited benchmark supernova models from sequence nlld2. (a) models nlld2r05tl0; (b) models nlld2rl0tl0; (c) models nlld2rl0tl5. 
Abrupt changes in the energy generation rate visible at t = 1 s (models nlld2r05tl0 and nlld2rl0tl0) and t = 0.75 s (models nlld2r05tl5) 
are caused by mesh derefinement. 



models, however these calculations were done in 3- 
dimensions. The energy generation rate in our models 
is initially smaller and increases at a faster rate than in 
the Garching models. The maximum generation rate is 
achieved at i w 1.25 s, roughly 0.5 s later than in the 
Garching models. In flash as well as in the Garching 
calculations the burning quenches ~ 0.4 s after the max- 
imum. 

Our model explosions are on average slightly more en- 
ergetic than those reported in Garching studies. The 
typical explosion en ergies are between pa 0.3 — 0.7 foe 
in ou r mod els while (iReinecke. Hillebran dt. fc Niemeveil 
120021 ) and (|R5pke fc Hillebrandt 2005a) reported explo- 
sions with energies « 0.35 and « 0.45, respectively. It 
is conceivable that the higher initial energy generation 
rates obtained in the Garching calculations may result 
in the lower explosion energies of their models since the 
faster initial expansion leaves less time for the flame to 
develop and burn stellar fuel. Such discrepancies can be 
explained by differences in both the adopted flame dy- 
namics model and the approximations used to describe 
the nuclear burning and not of serious concern in the 
context of the following results. 

3. DETONATING FAILED DEFLAGRATIONS 

Limited by the assumption of axial symmetry, we con- 
sidered off-center bubble ignitions in which the flame ini- 
tially occupies a small spherically symmetric region(s) 
positioned at the symmetry axis (r = cm). Two fami- 
lies of models were constructed, one with a single ignition 
point and the other with two ignition points. For the lat- 
ter, we consider only simultaneous ignitions, although in 
principle the multiple igni tion process could be e xtend ed 
in time (see, for example, ISchmidt fc Niemeveij (|2006D '). 

Following our verification study, off-center explosion 
models were obtained at the maximum resolution of 
8 km. In models with two ignition points, the flame 
regions were initialized in different hemispheres. Table [2] 
summarizes the parameters describing the initial flame 
and mesh configurations of the off-center ignition mod- 
els. Here Zb^i and Zf,^2 are the locations of the ignitions 
points along the symmetry axis (r — cm), and Rh 
is the radius of flame regions. To keep the flame re- 
gions well-separated, the mesh resolution in the central 



Model 



Y12 

Y25 

Y50 

YlOO 

Y70YM25 

Y100YM25 

Y75YM50 



TABLE 2 

Off-center Ignition Configurations 
zt,l [km] 2ft, 2 [km] Rj [km] 



12.5 

25 

50 

100 

70 

100 

75 



-25 
-25 
-50 



50 
50 
50 
50 
35 
50 
50 



300 km region was increased by a factor of 2 for a short 
period of time after ignition (0.4 s in models Y70Y1VI25 
and Y100YM25 and 0.1 s in model Y75YM50). This 
additional refinement was also necessary to adequately 
resolve the smaller bubbles {Rb — 35 km) used in model 
Y70YM25. 

The computational domain extended up to 131, 072 km 
and 524, 288 km in single- and double-ignition point mod- 
els, respectively. In anticipation of an extended and 
asymmetric evolution at early times, the region of adap- 
tive meshing was extended to 6,000 km in radius. The 
initial conditions did not include random velocity pertur- 
bations. All other simulation parameters were identical 
to those used in the verification study (see § 12.4. ip . 

3.1. Explosion phase 

than convenience and limits 
early multi-dimensional in- 
vestigations of white dwarf deflagrations assumed 
perhaps only slightly perturbed but otherwise spheri 
cally symmetric ignition conditions fMiiller fc Arnet 
1982t ,Livnc 1993; Arnctt & Livnc 1994a; Khokhlo_ 
19951 [20001: [Reinecke. Hillebrandt. fc Niemeveii 1200"^ 
Gamezo et al. 1 I2003D . Such a choice is not nec- 
essarily the most natural one given the white 
dwarf core is believed to be convective prior to 



For no other reason 
in computing power, 



11- 



runaway (iWooslevI fT990t iGarcia-Senz fc WooslevI 
[19951: iWooslev. Wunsch. fc KuhlenI 120041) . an ex- 
pectation supported by rec ent multi-dime n sional 
hydrodynamic investigations (iHoflich fc SteinI 120021 : 
'Kuhle n. Wooslev. fc Glatzmaieij|2006l ). This led several 
groups to consider progressively more complex and 
realistic (although not necessarily correct!) initial 
flame configurations ([Niemeyer. Hillebrandt. fc WooslevI 
IT99I IGarcia-Senz fc Bravd 120051 : iRopke et all 12006 : 
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ISchmidt fc Niemeveij[200i) . 

Here we adopt a similar approach, but not being dis- 
couraged by the failure of an initial deflagration to pro- 
duce a supernova, we continue our investigations through 
the following stages of evolution. Our preliminary inves- 
tigations (PCL) indicated that the energy released in the 
deflagration may be used to compress the stellar surface 
layers thereby forming seed points for detonations. We 
were unable, however, to study that process in detail 
at that time, and only speculated about the possibility. 
Here, we revisit our original idea of a deflagration to det- 
onation transition following a slightly off-center ignition. 

3.1.1. Failed deflagration phase 

In all models ignited off-center, the evolution initially 
proceeds in a way much similar to that described in 
PCL. Owing to strong buoyancy and the relative slow- 
ness of laminar burning, the whole burning region is 
quickly expelled from the core (iGarcia-Senz &: WooslevI 
ll995tlWooslev. Wunsch. fc Kuhlenll2004l ) consuming onlv 
a small amount of fuel on its way to the stellar surface. 
We describe the early evolution of burning regions in the 
set of single ignition models in terms of mean velocities, 
the position of centroids, and the effective radii of burned 
matter. The effective radii correspond to spheres of the 
same surface area as the burned region. 

From their original positions, the initially motionless 
bubbles are driven by buoyancy and develop primar- 
ily vertical velocities, as can be seen by comparing the 
mean vertical (Fig. [TOT a)) and radial (Fig. [rOT b)) veloc- 
ities. The initial acceleration is stronger in cases where 
the bubble is placed farther away from the stellar cen- 
ter (models Y50 and YlOO, thick and thin dashed lines in 
Fig.fTO]) and weaker when the bubble originates deeper in 
the core (models Y12 and Y25, thick and thin solid lines 
in Fig. nop . Although the mean motion of the burned 
region is away from the stellar center, the rich and com- 
plex flow field includes downflows and outflows develop- 
ing inside the region that lead to intermittent (apparent) 
deceleration and acceleration (seen as mild wiggles su- 
perimposed on the velocity curves). Each such wiggle in 
velocity is associated with the destruction of the current 
generation of Rayleigh- Taylor bubbles and the formation 
of the next. Our data indicates that perhaps two genera- 
tions of such bubbles are created during the deflagration 
phase. 

A sudden drop in vertical velocity and a rapid increase 
of lateral expansion marks the moment of bubble break- 
out. This phase is not well-defined but occurs roughly 
at « 0.7 s in model YlOO and not until t « 1 s in model 
Y12. The effective radius of the burned regions (defined 
as a radius of a sphere with surface area equal to the 
flame surface area) at breakout is very similar between 
the models, « 400 - 450 km (Fig. [TOfc)). However, the 
ccntroid of the burned region is located much farther out 
in model YlOO than in model Y12. That is understand- 
able considering there is more time for the RT instability 
to develop structure and for the region to grow laterally 
when the ignition takes place closer to the stellar center. 
This also has profound consequences for the evolution 
of the star. A longer deflagration phase allows for more 
burning, causes more matter being lifted from the stel- 
lar core, and eventually results in stronger stellar expan- 
sion during the early (t < 2 s) post-breakout evolution 



(Fig. [TTJ In the case of double point ignitions, the early 
evolution proceeds in a very similar way to that of single 
ignitions with a proportionally increase in the amount of 
burned material and stronger stellar expansion. 

3.1.2. Detonation phase 

One of the motivating factors behind extending our 
study to double ignition point scenarios was to examine 
whether a detonation can be formed when the maximal 
(and perhaps even boosted by flawed numerics) focusing 
offered by the symmetry axis is not present. Although 
we observed detonations forming in all off-center ignited 
models, only in one double ignition model, Y100YM25, 
does the detonation form near the equatorial plane. In 
the remaining two double ignition models, detonations 
eventually emerge near the symmetry axis. Although 
both models eventually detonate, we cannot consider 
them as examples of successful asymmetric DFDs. Nev- 
ertheless, both of them provide additional evidence for 
shock to detonation transition. In what follows we will 
first overview the formation of detonations in single ig- 
nition models. Then we will discuss the flow dynamics 
leading to detonation in model Y100YM25. We will con- 
clude by presenting the ejecta morphology soon after the 
shock breakout along stellar surface is complete. 

The progenitor structure around the time when a det- 
onation forms is shown in Fig.[T2]for models Y12, YlOO, 
and Y100YM25. In all models the bulk of progenitor has 
retained its original spherical characters. We do not find 
any substantial large-scale deformations, except that the 
stellar cores in single-ignition models are slightly ellip- 
soidal in shape with axis ratio 1.2 — 1.3. Stellar expansion 
has decreased the core density to ~ 9 x 10^ g cm~'^ in 
model Y12 and « 4 x 10^ g cm"^ YlOO. This is consis- 
tent with our expectation that lower central densities are 
to be found in models that experienced more energetic 
deflagrations. 

The stellar core is surrounded by an extended strongly 
turbulent atmosphere. Comparing three models shown 
in the upper row in Fig. I12i the atmosphere appears bet- 
ter developed (more extended and turbulent) in models 
that release more energy in the deflagration. This at- 
mosphere formed following the breakout of deflagration 
products through the stellar surface, at which time the 
ashes accelerated unburned surface layers both radially 
and laterally. This circular wave carried both fuel and 
products of the deflagration along the surface of the star. 
The following evolution depends on whether there was 
one or more ignition points. 

If the single ignition case, the surface wave eventually 
completely engulfs the progenitor and collides with itself 
in a region located opposite breakout. A conical shock 
wave forms in the process that thermalizes the kinetic 
energy of the incoming flow. This shock can be seen as 
a vertical structure near r = cm extending down from 
(r, z) « (8 X 10^ -3.5 x 10^) cm in model Y12 (Fig.IHlJa)) 
and (r, z) « (4 x 10^ -2.5 x 10^) cm in model YlOO 
(Fig.inrc)). 

In the multi-point ignition case, there will presumably 
be several breakout points and related surface waves that 
will be colliding with one another. Therefore, it is con- 
ceivable that several shock-dominated regions might be 
formed. Some of those shocks might be weaker and oth- 
ers stronger than the ones found here. The number of 
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Fig. 10. — Kinematics and growth of tlie burned region in a set of single ignition off-center deflagrations, (a) mean vertical velocity; (b) 
mean radial velocity; (c) equivalent radius; (d) vertical centroid position. Data for models Y12, Y25, Y50, and YlOO are shown with thick 
solid, thin solid, thick dotted and thin dotted lines, respectively. 



possible scenarios and outcomes is certainly much larger 
than represented in a limited sample of the initial con- 
figurations considered in this study. Nevertheless, we ex- 
pect that our models capture the essential features of the 
evolution. In the case of model Y100YM25, for example, 
the material of two waves collides near the equatorial 
plane forming a jet-like radial inflow and outflow near 
(r,z) « (2.5 X 10^,-5 x 10^) cm (Fig. [HJc)). This is 
essentially the same flow configuration we found in the 
single point ignition models. 

In what follows, we first focus our discussion on details 
of the transition to detonation process in three selected 
models. Then we characterize the evolution of the ex- 
ploding star during the passage of the detonation wave. 

Shock-to-detonation Transition — In all cases considered 
in our study, we found explosions foll owing a shock- 
to-detonation tran sition, SDT (see, e.g. JBdzil &: Kapilal 
119921 : ISharpeir2002l and references therein) . Although re- 
gions forming detonations differ greatly in morphology. 



the common ingredients of the process include the pres- 
ence of a strong acoustic wave, dense fuel, and a pro- 
longed compression of the region. This is illustrated in 
the bottom row of Fig. [T^ which shows the density distri- 
bution and major flow structures involved in a transition 
to detonation process. 

Model Y12 shares the initial conditions with the origi- 
nal PCL study. In the Y12 model, we found no sign of a 
possible transition to detonation for t < 3.5 s (Fig. fl^ a)) 
in contrast to the PCL study, in which a detonation wave 
formed at i « 1.9 s This difference in timing is due solely 
to the incorrect energetics of approximate deflagration 
burner used in the PCL calculations. The overestimate 
of the energy release and hence buoyancy by a factor of 
« 3 in the PCL model resulted in a much shorter de- 
flagration phase, a lower overall energy release, a more 
compact progenitor and ultimately a significantly ear- 
lier formation of the detonation. At the same time, the 
less expanded progenitor allowed for the surface wave to 
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Fig. 11. — Evolution of the effective stellar radius (defined as 
the radius of a sphere with the same volume occupied by matter 
with density > 1 X 10* g cm~^) in the single ignition off-center 
deflagrations. Data for models Y12, Y25, Y50, and YlOO are shown 
with thick solid, thin solid, thick dotted and thin dotted lines, 
respectively. 

move at a relatively higher speed (due to a lower orbit) 
and with higher post-shock temperatures, possibly en- 
hancing the likelihood of a transition to detonation. 

Similar to PCL, we observe the formation of a conical 
shock in Y12. However, the conditions inside the shocked 
region allow only for residual burning that perturbs the 
shocked gas (the low density region near symmetry axis 
at r = cm in the lower section of Fig. [T^ a)). The 
wave that transits to detonation is an accretion shock 
(blue isocontour extending horizontally from (r, z) « 
(0, -2.55 X 10^) cm in Fig. [I2i;d)) created by infalhng 
material which is trapped horizontally by the symmetry 
axis and the incoming deflagration products, and verti- 
cally by the bulk of the stellar material and the stagna- 
tion point formed behind the conical shock. A transition 
to detonation takes place at the symmetry axis where 
the ram pressure of the incoming flow and the resulting 
post-shock temperature are the highest. (Although the 
accretion flow is predominantly along z-direction, there 
exists a lateral velocity gradient, AvrjAr ~ —3, in the 
flow that focuses the flow toward the symmetry axis.) 

In model YlOO (Fig. [T^ d)). a detonation wave can 
be seen as a nearly vertically propagating shock located 
near (r, z) « (6.5 x 10'', —2.85 x 10^) cm. Unlike in model 
Y12, here the detonation is not directly associated with 
the symmetry axis. The conical shock is visible as an 
wave originating near (r, z) « (7.5 x 10*", —3 x 10*) cm, 
just to the right of the detonation region. The appar- 
ent closeness of the flame front is coincidental and has 
nothing to do with origin of the detonation wave. By the 
time of Fig. fl^ d) the conical shock is already sweeping 
through deflagration products, however the detonation 
appears moving through a channel of shocked fuel even- 
tually connecting to the bulk of stellar material. 

Model Y100YM25 displays by far the most complex 
flow structure in the region where a transition to deto- 
nation occurs. The configuration of two colliding surface 



waves resembles that of a "self-colliding" surface wave 
found in models with a single ignition point. Here, how- 
ever, a symmetry of the problem is broken. Not only 
does the collision not take place at the symmetry axis, 
but the ignition points were initially located at differ- 
ent distances from the core. The timing, energetics, and 
morphology of each wave were therefore slightly differ- 
ent. This difference eventually leads to a shift of the 
collision plane « —400 km from the equator (Fig. [T2r f)). 
Furthermore, the collision does not occur "head-on" but 
rather material from the slower wave (ignition point lo- 
cated closer to the core) tends to penetrate underneath 
that of the faster wave (ignition point located farther 
away from the core). In the end, the whole region shows 
a slight tendency to roll. 

Another difference from the highly symmetric single 
ignition models is that the broken symmetry offers the 
potential for creating more than just one shocked re- 
gion. This is indeed the case in model Y100YM25. 
Two shock fronts moving in the radial direction can 
be seen inside the collision region: one located near 
(r, z) w (2.3 x 10^,-4 x 10*") cm, and another near 
(r, z) « (2.1 X 10^,-4 X 10'') cm. Both fronts are cre- 
ated near the collision plane which might be understood 
given this is where we expect the thermalization rate 
of the colliding flows being the greatest. The former 
shock wave evolves into a self-sustained detonation while 
the latter soon dies off. Once a detonation is formed, 
the wave travels approximately along a fuel-rich channel 
(a flame-bounded horizontally extending structure near 
(r, z) « (2.5 X 10^,-3.5 x lO"") cm that connects to the 
bulk of the unburned stellar material. 

Some more details and observations can be offered re- 
garding the SDTs observed in our subset of models. We 
found that transitions to detonation occur in gas with 
pre-shock densities w 1 — 3 x 10^ g cm^'^ in models Y12 
and YlOO and > 5 x 10^ g cm'^ in model Y 100YM25. As 
demonstrated by (jArnett fc Livnelll994bf ). at these den- 
sities the typical radius of region that can successfully 
transit to detonation is '^few kilometers. This is smaller 
than the numerical resolution in our models. It is con- 
ceivable that problems with producing SDTs in some of 
our double-point ignition models might be attributed to 
insufficient resolution. For the same reason, the observed 
SDTs may require less time to launch detonations after 
a strong shock wave forms. However, preconditioning 
of the fuel for SDT is certainly a temporally extended 
process requiring both compression of the material and 
thermalization of the flow. The latter is aided by the 
confinement that makes the thermalization process more 
efficient. Still, large amounts of energy need to be sup- 
plied to the region and the typical velocity jumps across 
shock waves are 4, 000 — 6, 000 km s~^. This guarantees 
post-shock temperatures > 1 x 10^ K, high enough for 
nuclear burning to take control of the flow dynamics. 

Evolution through detonation — Transitions to detona- 
tions occur in different models at different times and 
locations although, as we discussed earlier, several nec- 
essary elements (high density fuel-rich matter, strong 
wave, extended confinement of the region) are commonly 
present. Figure [T2] shows the evolution of the equiva- 
lent stellar radius (defined as the radius of a sphere with 
the same volume as that occupied by matter of density 
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Fig. 12. — Post-deflagration progenitor structure around the time of transition to detonation. Shown is density in log scale, (left): model 
Y12, (a) t = 3.575 s, (b) t = 3.6 s; (middle): model YlOO, (c) t = 2.0 s, (d) t = 2.125 s; (right): model Y100YM25, (e) t = 1.776 s, (f) 
t = 1.8 s. Several structures are shown with isocontours in the lower row: T = 1 X 10® K (red), flame front (green), shocks (blue). 



> 1 X 10'* g cm~'^) in our sample of DFD models. This 
initial steady stellar expansion is due to energy deposi- 
tion by a failed deflagration. In most cases a detonation 
occurs when the progenitor either approaches or begins 
to recollapse. This is understandable since at later times 
the energy of surface waves is likely to quickly dissipate, 
thus decreasing the likelihood of strong hydrodynamic 
interactions taking place (e.g. in multi-point ignitions, 
model Y100YM25). Alternatively, the accretion flows 
can develop only once expansion stops and high accre- 
tion luminosities (required for SDT) can be expected 
only shortly after accretion flows develop. In particular, 
in models Y70YM25 and Y75YM50 (thick dotted and 
dotted hue in Fig. [131 respectively) no SDT occurred 
during the collision of surface waves. Instead, flow per- 
turbations accumulated in the regions near the symmetry 
axis, evolved into jet-like flows and eventually triggered 
detonations. With our verification study providing evi- 
dence that the evolution of perturbation near the sym- 
metry axis cannot be entirely trusted, we do not consider 
these two models as successful DFDs produced by dou- 
ble point ignitions. (One probably could still consider 
them members of single point ignition families, perhaps 
obtained from different initial conditions, and provided 
their radii prior to SDT were similar to original single 
ignition models. We do not consider this inelegant pos- 



sibility any further.) 

Figure [T3] shows the morphology of the exploding 
model supernovae Y12, YlOO, and Y10GYM25 shortly 
after the central density drops below 1 x 10® g cm~'^ and 
nuclear burning essentially quenches. Several important 
observations can be made. All model supernova ejecta 
are stratified. The ejecta are composed of a featureless 
core surrounded by inhomogeneous external layers. This 
composite structure of the ejecta is a direct result of the 
two stages of the explosion, each involving a diametri- 
cally different dominant process. The core is the rem- 
nant of the original progenitor which has been expanded 
by energy released during the defiagration phase, the 
density distribution in the core displays slight asymme- 
try reflecting the character of the initial conditions (the 
shift in the density maximum to the lower hemisphere 
present in single ignition models is absent in double ig- 
nition model) . The outer layers may show global asym- 
metry, especially in models with a single ignition point, 
but are always rich in structures reflecting the turbulent 
nature of the deflagration and the violent evolution of 
surface waves. 

The density stratification of the ejecta is accompanied 
by compositional stratification. This is due to the deto- 
nation wave which synthesized iron peak elements when 
sweeping through the dense central regions and interme- 
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Fig. 13. — Evolution of effective stellar radius (defined as ra- 
dius of a sphere with the volume occupied by matter with den- 
sity > 1 X 10* g cm~^) in model Y12 (thick solid), Y25 (medium 
solid), Y50 (solid), YlOO (thin solid), Y70YM25 (thick dotted), 
Y100YM25 (medium dotted), and Y75YM50 (dotted). Evolution 
toward SDT proceeds through overall continuous expansion fol- 
lowed by a possible period of recoUapse. Shock breakout is followed 
by a very rapid expansion of stellar material at approximately con- 
stant velocity. 

diate group elements when encountering the less dense 
outer layers. The detonation is driven by essentially in- 
stantaneous energy release due to carbon burning fol- 
lowed by approach to NSQE (silicon-group elements) and 
final relaxation to NSE (iron-group elements) . With rel- 
atively crude resolution, only the relaxation to NSE and 
approach to NSQE at the lowest densities c an be con- 
sider ed spatially resolved in our simulations (|Khokhlovl 
119891 Fig. 9). This problem, however, is not related to 
the energy release and so does not influence the overall 
dynamics of the detonation front. The low resolution of 
our models prevents us from considering possible effects 
related to the curvature of the detonation front. Once 
again, such effects affect the s tructure of th e detonation 
wave only on scales < 10 km ()Sharpdl2001h and are not 
expected to influence the large scale dynamics of the det- 
onation wave. 

The outermost ejecta are composed of unburned fuel 
mixed with the deflagration products, most likely inter- 
mediate group elements with locally entrained iron-rich 
material. Close to the core, those layers were overrun 
by the detonation wave that additionally modified their 
composition. Although definitely present in our calcula- 
tions, we estimate this effect to be small. To calculate the 
detailed composition of the ejecta, including the deflagra- 
tion miX j requires p ostproce ssing nucl eosynthesis (see for 
example iTravaglig Hillcbra ndt. fc Re inecke 2006) and is 
beyond the scope of this presentation. 

Detonations are known to be susceptible to trans- 
yerse perturbations a nd developing cellular structure 
(|Fickett fc Davis! 120011 Chap. 7). Potential sources for 
such perturbations are abundant in our models and in- 
clude upstream flow perturbations due to preexisting 
convection and turbulence (not considered here), initial 
deflagration, or possible numerical oscillations of grid- 



aligned shock fronts, i.e. odd-even decoupling (jQuirkl 
I1994D . We found no clear evidence for cellular structure 
in our calculations. One possibility is that the numeri- 
cal resolution i s insufficient to resolve detona tion cells of 
sub-km size ( Gamezo et al.lll999l : ll^ [200(7). Also, the 
time available for cellular structure is very limited as low 
density material rapidly expands following the passage of 
the detonation wave and the burning quickly quenches. 

The evolution of the explosion energy (kinetic -I- in- 
ternal -|- potential) and burned mass from the moment 
of ignition until the shock breakout is shown in Fig. 1151 
With plenty of unburned fuel available to the detonation, 
DFD supernovae are energetic events with typical post- 
detonation energies « 1.3 — 1.5 foe. The deflagration 
phase typically supplies only 0.06 to 0.15 foe of energy 
in burning < 0.1 Mq of stellar fuel. The bulk of the en- 
ergy is released during about 0.4 s when the detonation 
wave sweeps through the white d warf. These findings re- 
semble the results presented by (jArnett fc Livneiri994bL 
Fig. 3) although the mechanism behind the transition to 
detonation is different in the two models. 

3.2. Homologous expansion: final properties 

The direct comparison of hydrodynamic explosion 
models to observations is accomplished through the cal- 
culation of synthetic model light curves, spectra, and 
possibly also spectrum polarization. These radiative 
transfer calculations take as an input the model super- 
nova ejecta with complete specification of the density 
distribution, ejecta chemical composition, and typically 
make the simplifying assumption that velocity linearly 
increases with distance from the ejecta center. This 
last assumption is satisfied to different degrees in vari- 
ous parts of the ejecta and generally does not hold true 
during the early stages of supernova expansion. The rea- 
son is that linear expansion requires establishing a fine 
balance of accelerations between neighboring fluid ele- 
ments in the ejecta, and this requires time. For example, 
our post-detonation explosion models contain large re- 
gions where energy of the flow is dominated by internal 
energy. This indicates the potential of fluid elements 
performing some work, possibly adjusting their motion 
relative to their neighbors. To allow for that process to 
operate and establish homologous expansion, the post- 
explosion needs to be continued for an extended period 
of time. Detailed discussion of this process in the context 
of Typ e la su pernova modeling was recently presented by 
iRopkel ([2005h . 

We obtained a complete set of homologously expanding 
model ejecta using flash and its adaptive mesh capabil- 
ities. Post-detonation models were carefully transported 
to a high-resolution uniform mesh with typical relative 
errors of total mass, total energy, and abundances not 
exceeding 0.1%, 0.5%, 5%, respectively. Given the sev- 
eral sources of uncertainties and variations in the original 
models, this accuracy is sufficient for any practical pur- 
poses. The interpolated models were then used to define 
the initial conditions in the flash calculations. 

We used a ratio of kinetic energy to the sum of inter- 
nal plus gravitational energies to monitor the approach 
of ejecta to homology. In course of several trial compu- 
tations, we established that continuing calculations for 
100 s after explosion guarantees that the energy ratio > 
100 anywhere in the ejecta except for the innermost half 
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Fig. 14. — Post-detonation structure of the exploding supernova at the time when the central density drops to 1 X 10^ g cm ^ and 
burning effectively quenches. Panels (a)-(c) show the density in log scale in models Y12 (i = 4.90 s), YlOO {t = 3.375 s), and Y100YM25 
(i = 2.878 s). Abundance isocontours X(^®Ni) = 0.95 and X(^®Si) = 0.2 are shown with black and gray lines, respectively. 
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Fig. 15. — Evolution of the total (explosion) energy through the 
end of detonation phase in model Y12 (thick solid), Y25 (medium 
soUd), Y50 (solid), YlOO (thin soUd), Y70YM25 (thick dotted), 
Y100YM25 (medium dotted), and Y75YM50 (dotted). Only small 
amount of energy, 0.06 — 0.15 foe is released during deflagration 
phase. Detonation phase lasts < 0.4 s. 

of the central nickel-rich core (and parts of essentially 
massless shocked ambient medium) . This required using 
a computational grid extending to 1.68 x 10^^ cm. Cal- 
culations were performed using an automatic mesh dere- 
finement scheme that kept the computational cost ap- 
proximately constant as the ejecta expanded. The effec- 
tive mesh size varied from 8, 192 initially to 2, 048 — 4, 096 
zones per dimension at a final time. 

The density and compositional structure of the homol- 
ogously expanding model ejecta are shown in Fig. 1161 
Here we show only the innermost part of the grid with 
expansion velocities reaching « 50, 000 km s""'^ at r = 
5 X 10^^ cm. The remaining part of the volume contains 
a low-density shocked ambient medium and the super- 
nova shock. The bulk of the ejecta material displays es- 
sentially the same structure as in early post-detonation 
models discussed in the previous section. Most visi- 
ble differences can be found in the outermost regions 



that here were already swept by the shock. The com- 
posite character of the bulk of the ejecta is preserved 
with a featureless core rich in iron group elements and 
the outer strongly mixed layers rich in silicon group el- 
ements. The compositional dichotomy of the outer lay- 
ers also reflects the contribution of two processes to the 
explosion. The inner well-deflned silicon-rich ring also 
contains an admixture of calcium (shown with white iso- 
contours in Fig.[Tni). This is due to the detonation wave 
nucleosynthesis calculation done with aproxl3t nuclear 
network. The outer silicon-rich and rather fragmented 
shell is devoid of calcium as this species was not consid- 
ered in the approximate 3-stage flame burner. Improving 
upon the approximate nucleosynthesis is one of the ur- 
gent future tasks, especially given the importance of the 
outer ejecta layers in formation of supernova spectrum 
(|Kasen fc Plewal[2005[ ) . It is also interesting to note that 
the each deflagration region seems responsible for form- 
ing its own outer silicon-rich ring. This is evident in 
model Y100YM25 (Fig. [HJc)). 

Table [3] presents approximate nucleosynthetic yields 
and flnal explosion energies for the complete set of ho- 
mologous DFD models. The homologous character of the 
models is confirmed by the consistently small fraction of 
potential and internal energies as compared to the total 
energy. In most models, explosion energies are in the 
range 1.3 — 1.5 foe. These produce between 0.9 to al- 
most 1.2 solar masses of ^^Ni and between 0.1 and 0.24 
solar masses of intermediate elements. Although model 
nickel masses may appear relatively high at first, such 
high nickel mass es might be typical for significant frac- 
tion of objects (jStritzinger et al.1l2006D . Furthermore, 
our estimates of nickel mass are likely the upper lim- 
its given aproxl3t nucleosynthesis does not account for 
production of other iron-group elements, e.g. stable iso- 
topes like ^"'Fe. Model Y75YM50 is the least energetic 
{Et « 1.08 X 10^^ erg), produces the least amount of 
nickel (« 0.51 Mq) and more than a half solar mass of 
intermediate mass elements. As we discussed earlier, we 
believe that this model should not be considered as a 
DFD, as the shock to detonation transition was likely 
being promoted by numerics. 

4. DISCUSSION 
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Fig. 16. — Homologously expanding model ejecta ^ 100 s after explosion. Panels (a)-(c) show density in log scale in models Y12, YlOO, 
and Y100YM25. Abundance isocontours X(56Ni) = 0.95, X(2**Si) = 0.2, and XC^^Ca) = 0.05 are shown with black, gray, and white 
lines, respectively. Expansion velocity is f» 50,000 km s~^ at distance 5 X 10^^ cm from the ejecta center. Only the innermost part of 
computational grid is shown. 



TABLE 3 
Homologous DFD models^ 



Model 

Ei 

—Ep 

4He 

12 Q 

16 o 

20Ne 

24Mg 

28 Si 
32s 

36 Ar 
«Ca 

44^- 

48 Or 
52pg 

56Ni 



Y12 



Y25 



Y50 



YlOO 



Y75YM25 Y100YM25 Y75YM50 



1.357 
1.59X 
2.52X 
8.03X 
8.73X 
0.107 
4.41 X 
8.70X 
0.127 
7.03X 
1.64X 
1.82X 
1.41X 
2.96X 
6.50X 
0.926 



10-* 
10-3 
10-3 
10-3 



1.496 

8.38x10-5 
2.39x10^3 
1.13x10-2 
5.49x10-3 
4.65x10-2 
3.79x10-'' 
3.40x10-2 
7.28x10-2 
3.65x10-2 
8.26x10-3 
8.95x10-3 
9.35x10-6 
1.49x10-"' 
3.43x10-3 
1.147 



1.515 

7.15x10-5 
2.38x10-3 
1.15x10-2 
3.30x10-3 
4.48x10-2 
3.28x10-'' 
3.42x10-2 
6.07x10-2 
3.06x10-2 
6.91x10-3 
7.53x10-3 
3.02x10-5 
1.42x10-" 
3.01x10-3 
1.173 



1.516 

7.09x10-5 
2.38x10-3 
1.10x10-2 
4.56x10-3 
3.91x10-2 
4.78x10-* 
2.81x10-2 
5.74x10-2 
3.18x10-2 
7.36x10-3 
8.09x10-3 
1.35x10-5 
1.42x10-* 
2.91x10-3 
1.186 



1.464 

5.34x10-* 
2.31x10-3 
1.03x10-2 
1.29x10-2 
7.54x10-2 
1.04x10-3 
4.51x10-2 
8.00x10-2 
4.21x10-2 
3.97x10-3 
1.02x10-2 
2.71x10-5 
1.78x10-* 
3.49x10-3 
1.075 



1.384 
2.87x10-5 
2.30x10-3 
8.36x10-3 
2.05x10-2 
9.82x10-2 
9.53x10-* 
6.74x10-2 
0.137 

8.75x10-2 
2.07x10-2 
2.20x10-2 
2.71x10-5 
3.43x10-* 
6.85x10-3 
0.895 



1.075 

1.97x10-3 
2.56x10-3 
2.25x10-3 
2.52x10-2 
0.237 

9.73x10-* 
0.194 
0.202 
0.124 

2.95x10-2 
3.24x10-2 
2.58x10-5 
4.83x10-* 
1.03x10-2 
0.510 



We studied the fate of a massive carbon/oxygen white 
dwarf foUowing an off-center mild ignition. We found 
that such initial configurations do not produce direct ex- 
plosions. Only a small amount of stellar fuel is initially 
consumed and the released energy is used to expand the 
progenitor. This is in agreement with several previous in- 
dependent studie s in which th e deflagration was either in- 
trinsically weak (|Arnett fc Livne 1994bi) or was initiated 
slightly off-center (iNiemever. Hillebrandt. &: WooslevI 
fl996t iLivne. Asida. fc Hoflichl l2005^. 

We found that the following evolution of the perturbed 
stellar material leads to the formation of isolated wave- 
dominated regions inside unburned material. We consid- 
ered these regions capable of launching a detonation wave 
through a shock-to-detonation transition. We observed 
the resulting detonations eventually consuming the bulk 
of the unburned progenitor. These detonating failed de- 
flagrations are energetic events with explosion energies 
« 1.3 - 1.5 foe. 

The model DFD ejecta appear composite, reflecting 
the presence of two different physical processes contribut- 
ing to the explosion. The central parts of the ejecta are 
composed of a mildly deformed but completely feature- 
less central core rich in iron peak elements. Stronger 
deformations may require different physics, e.g. rotation 



(jHoflichll2005[) . The core is surrounded by a slightly in- 
homogeneous inner ring rich in silicon group elements, a 
product of the detonation burning at low densities. Fi- 
nally, the outermost layer is highly turbulent containing 
a mix of deflagration products and unburned material. 
Preliminary nucleosynthesis results indicate that DFD 
models typically produce over 0.9 Mq of iron group el- 
ements and 0.3Mq intermediate elements. The burn is 
almost complete leaving essentially no carbon. 

Our conclusions are based on calculations using a re- 
vised numerical scheme that contains substantial im- 
provements. We found that the energet ics of defiagra- 
tion stage originally considered in (C alder et al.l I2004L 
PCL) tended to overestimate buoyancy effects by a fac- 
tor « 3. We used a set of self- heating nuclear network 
calculations and implemented a density-dependent en- 
ergy release scheme. Additional modifications to the 
approximate nucleosynthesis were included to improve 
the dynamics of the early phases still further. With 
the revised energetics, we calibrated the numerical flame 
speed to match the resul ts of detailed calculations of 
iTimmes fc WooslevI (|1992f ). 

The revised deflagration code was subsequently ver- 
ified against an independent set of results of centrally 
ignited deflagrations obtained with the Garching super- 



Dynamics of Detonating Failed Deflagrations 



19 



nova code. Good qualitative agreement was found be- 
tween the two codes. The database of computer models 
is offered on-line to facilitate future verification (code-to- 
code comparison) studies. 

Furthermore, a numerical procedure to stabilize model 
progenitors has been developed. These stabilized progen- 
itors not only provided initial conditions for supernova 
simulations but also allowed us to examine the fidelity of 
hydrodynamic advection in axisymmetric situations. In 
particular we found that, on the one hand, no numeri- 
cally stable progenitors can be obtained if resolution is 
too low and, on the other hand, small perturbations are 
strongly amplified near the symmetry axis in highly re- 
solved models. This analysis allowed us to identify the 
optimal resolution for our supernova calculations. 

We analyzed the observed pattern of shock-to- 
detonation transitions (SDTs) in some detail. We iden- 
tified the presence of sufficiently dense fuel, strongly ki- 
netic flow, wave formation, and a persistent confinement 
of the region with additional pressure increase due to nu- 
clear burning in the shocked gas as necessary conditions 
for a SDT. We found that the phenomenon of SDT is 
not exclusively associated with the presence of a symme- 
try axis. We also found that SDTs can occur in regions 
completely free of possible geometrical boosting effects, 
i.e. near the equatorial plane. 

However, transitions to detonations were not a robust 
prediction of such models. There are some possible rea- 
sons for that. For example, by assuming axisymmetry we 
eliminate an angular direction in which additional per- 
turbations may develop. Such perturbations will enrich 
the flow field creating more seed points for SDT and, 
at the same time, increase the amplitude of fluctuations. 
That is, the assumed symmetry is likely limiting the pos- 
sible wave interactions and presumably denying extreme 
events such as a SDT. In addition, for stability reasons, 
our calculations had to be performed using suboptimal 
numerical resolution. This caused a strong numerical 
damping and limited sampling of the perturbed regions 
harming chances for observing SDTs still further. On 
the other hand, we did not include physics that may, ef- 
fectively, make the system behavior appear more viscous 
(e.g. magnetic flelds) inhibiting formation of small scale 
structures. 

Our findings also hinge on the assumption that the 
SDT process is relatively insensitive to details of evolu- 
tion on scales unresolved in our simulations. This re- 
mains to be demonstrated, ideally in the course of dedi- 
cated highly-resolved model calculations of compression- 
ally heated fuel-rich degenerate mixtures. It will be 
a daunting task. If any parallel can be drawn, expe- 
rience accumulated by modelers of incrtially confined 
fusion systems might be of great help in su ch studies 
(jAtzeni fc Mever-ter-Veiml 120041 : lOrakg l2006f ) . Even if 
such models are successfully constructed, many doubts 
will remain regarding the outcome of such calculations 
given how limited our knowledge about real systems is. 
For example, as we mentioned before (PCL), although we 
consider a pure carbon/oxygen progenitor it is almost 
certain that in nature p rogenitor's surface layers con- 
tain ad mixture of helium (|Nomotalll982t lYoon fc Langed 
l2005b( l. Compositional changes will affect energetics of 
nuclear burn adding entirely new dimension to the prob- 
lem. 



Being mindful of numerous simplifying assumptions 
and model inaccuracies, the essential findings of this 
work are, therefore, rather modest and can be sum- 
marized as evidence of strong, localized, and prolonged 
shock heating in regions rich in fuel. We note that these 
are necessary conditions for shock-to-detonation transi- 
tion to occur. We believe this observation is independent 
of particular details of our model, especially numerics, 
making SDT one of prime suspects for triggering deto- 
nations in SN la. 

However, even if no seed point forms a detonation 
through SDT, this second, after the initial deflagration, 
failure to unbind the star in no way automatically im- 
plies supernova will not occur. Perhaps just the opposite. 
The extensive large-scale mixing of deflagration products 
with unburned outer stellar layers combined with abun- 
dantly present strong acoustic perturbations appear the 
conditions are ripe for the Zeldovich gradient mechanism 
(jKhokhlov. Oran. fc Whed^ [l99l. We_ may ex pect 
that for rotating progenitors ijYoon fc Langeij2005a| ) per- 
turbations of surface layers will be even stronger due to 
presence of additional shear component. If all these op- 
portunities are missed, the white dwarf might still be 
given another chance to produce the supernova. A failed 
attempt to explode would then be a beginning of a cy- 
cle that repeats, perhaps several times, as the expanded 
white dwarf eventually cools down, shrinks, and prepares 
for another ignition. That is, the explosion process might 
be a lengthy one, a kind of laborious slow death. 

As the observations improve, we are also beginning 
to collect evidence that the degree of diversity of SN la 
might be greater than original anticipated (or desired!). 
Re cent observations of the pecuhar supernova SN 2002cx 
by I.Tha et all (120061 ) argue in favor of low energy ex- 
plosion and large degree of mixing in this object, two 
characteristics that essentially preclude any involvement 
of a det onation in the ex plosion process. Other objects 
hsted bv lJha et al] (|2006f ) may belong to SN 2002cx class. 
These rare peculiar supernovae might be genuine exam- 
ples of pure deflagrations. Or perhaps these are objects 
that underwent several failed deflagrations leaving only 
small amount of material to fuel the detonation. If so, 
normally bright supernovae might be DFDs that suc- 
ceeded early, with the occurrence of a detonation (or lack 
thereof) being the primary element determining the ob- 
servational properties of a given event. 

On the other hand, SNe la display some characteristics 
that we find difficult to explain in the framework of pure 
defiagrations. O ne example are iron-rich clumps found 
in Ty pe la SNRs (jWarren fc Hughesll200-A IWarren et all 
l2005f l. In the model proposed here, the inner ring of 
intermediate elements seems to be a natural site for 
producing nickel -rich blobs that may float away from 
the central core (iBlondin. Borkowski. fc Revnoldj 120011 : 
IWang fc Chevaliei1l2001h . Those radioactive blobs may 
also be explained in pure deflagration models that nat- 
urally produce clumpy ejecta. However, dominant and 
isolated regions rich in iron-group elements like the one 
observed in Tycho SNR (IVancura. Gorenstein. fc Hughes! 
[19951 Irniren et all iMsF can hardly be produced in a 
pure deflagration in which several such regions are ex- 
pected to be simultaneously present. In DFDs, such an 
isolated clump located near the outer edge of the super- 
nova remnant might be a material burned deep in the 
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progenitor core and transported to the stellar surface by 
one of the deflagrating bubbles. 

We cannot address the above question without detailed 
nucleosynthesis calculations. This is one of the possi- 
ble future directions for research. In addition, the re- 
laxation of the assumption of axial symmetry, although 
costly, will be necessary. But even with only the current 
approximate nucleosynthesis and simplifled geometry of 
the problem, we are in a position to conduct prelimi- 
nary validation studies against observations for a subset 
of DFD models. This will be the subject of the next 
communication in this series. 
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